
Из всех приложений машинного обучения диагностика любого серьезного заболевания с помощью черного ящика всегда будет сложной задачей. Если результатом модели является конкретный курс лечения (потенциально с побочными эффектами), операция или отсутствие лечения, люди захотят знать почему.
Этот набор данных дает ряд переменных вместе с целевым состоянием наличия или отсутствия сердечного заболевания. Ниже данные сначала используются в простой модели случайного леса, а затем модель исследуется с использованием инструментов и методов объяснимости машинного обучения.
Импортируем соответствующие библиотеки:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns #for plotting
from sklearn.ensemble import RandomForestClassifier #for the model
from sklearn.tree import DecisionTreeClassifier
from sklearn.tree import export_graphviz #plot tree
from sklearn.metrics import roc_curve, auc #for model evaluation
from sklearn.metrics import classification_report #for model evaluation
from sklearn.metrics import confusion_matrix #for model evaluation
from sklearn.model_selection import train_test_split #for data splitting
import eli5 #for purmutation importance
from eli5.sklearn import PermutationImportance
import shap #for SHAP values
from pdpbox import pdp, info_plots #for partial plots
np.random.seed(123) #ensure reproducibility
pd.options.mode.chained_assignment = None #hide any pandas warnings
Загрузим данные Next, load the data,
dt = pd.read_csv("../input/heart.csv")
Выведем на экран первые 10 строк
dt.head(10)
Это чистый, простой для понимания набор данных. Однако значение некоторых заголовков столбцов неочевидно. Вот что они имеют в виду:
Чтобы избежать HARKING (или выдвижения гипотез после того, как известны результаты), взглянем на онлайн-руководства по диагностике сердечных заболеваний, и найдём некоторые из приведенных ниже терминов.
Diagnosis: Диагноз сердечных заболеваний ставится на основании сочетания клинических признаков и результатов анализов. Типы тестов будут выбраны на основе того, что думает врач 1, от электрокардиограммы и компьютерной томографии сердца (КТ) до анализов крови и нагрузочных тестов 2.
Просмотр информации о факторах риска сердечных заболеваний привел меня к следующему: высокий уровень холестерина, высокое кровяное давление, диабет, вес, семейный анамнез и курение 3. Согласно другому источнику 4, основные факторы, которые могут не подлежат изменению: возраст, мужской пол и наследственность. Обратите внимание, что талассемия, одна из переменных в этом наборе данных, - это наследственность. Основными факторами, которые можно изменить, являются: Курение, высокий уровень холестерина, высокое кровяное давление, отсутствие физической активности, избыточный вес и диабет. К другим факторам относятся стресс, алкоголь и неправильное питание.
Я не могу найти ссылки на статьи о «количество крупных сосудов», но, учитывая, что болезнь сердца определяется как «... что происходит, когда кровоснабжение вашего сердца блокируется или прерывается из-за накопления жирных веществ в коронарные артерии », кажется логичным, что больше крупных сосудов - это хорошо и, следовательно, снижает вероятность сердечных заболеваний.
Учитывая вышеизложенное, я бы предположил, что, если модель обладает некоторой способностью к прогнозированию, мы увидим эти факторы как наиболее важные.
Давайте изменим имена столбцов, чтобы они были немного понятнее:
dt.columns = ['age', 'sex', 'chest_pain_type', 'resting_blood_pressure', 'cholesterol', 'fasting_blood_sugar', 'rest_ecg', 'max_heart_rate_achieved',
'exercise_induced_angina', 'st_depression', 'st_slope', 'num_major_vessels', 'thalassemia', 'target']
Также изменим значения категориальных переменных, чтобы позже удобнее было объяснять:
dt['sex'][dt['sex'] == 0] = 'female'
dt['sex'][dt['sex'] == 1] = 'male'
dt['chest_pain_type'][dt['chest_pain_type'] == 1] = 'typical angina'
dt['chest_pain_type'][dt['chest_pain_type'] == 2] = 'atypical angina'
dt['chest_pain_type'][dt['chest_pain_type'] == 3] = 'non-anginal pain'
dt['chest_pain_type'][dt['chest_pain_type'] == 4] = 'asymptomatic'
dt['fasting_blood_sugar'][dt['fasting_blood_sugar'] == 0] = 'lower than 120mg/ml'
dt['fasting_blood_sugar'][dt['fasting_blood_sugar'] == 1] = 'greater than 120mg/ml'
dt['rest_ecg'][dt['rest_ecg'] == 0] = 'normal'
dt['rest_ecg'][dt['rest_ecg'] == 1] = 'ST-T wave abnormality'
dt['rest_ecg'][dt['rest_ecg'] == 2] = 'left ventricular hypertrophy'
dt['exercise_induced_angina'][dt['exercise_induced_angina'] == 0] = 'no'
dt['exercise_induced_angina'][dt['exercise_induced_angina'] == 1] = 'yes'
dt['st_slope'][dt['st_slope'] == 1] = 'upsloping'
dt['st_slope'][dt['st_slope'] == 2] = 'flat'
dt['st_slope'][dt['st_slope'] == 3] = 'downsloping'
dt['thalassemia'][dt['thalassemia'] == 1] = 'normal'
dt['thalassemia'][dt['thalassemia'] == 2] = 'fixed defect'
dt['thalassemia'][dt['thalassemia'] == 3] = 'reversable defect'
Проверимм типы данных:
dt.dtypes
Некоторые из них не совсем правильные. Приведенный ниже код преобразует их в категориальные переменные:
dt['sex'] = dt['sex'].astype('object')
dt['chest_pain_type'] = dt['chest_pain_type'].astype('object')
dt['fasting_blood_sugar'] = dt['fasting_blood_sugar'].astype('object')
dt['rest_ecg'] = dt['rest_ecg'].astype('object')
dt['exercise_induced_angina'] = dt['exercise_induced_angina'].astype('object')
dt['st_slope'] = dt['st_slope'].astype('object')
dt['thalassemia'] = dt['thalassemia'].astype('object')
dt.dtypes
Для категориальных переменных нам нужно создать фиктивные переменные. Я также собираюсь исключить первую категорию каждого из них. Например, вместо «мужчина» и «женщина» у нас будет «мужчина» со значениями 0 или 1 (1 - мужчина, а 0 - женщина).
dt = pd.get_dummies(dt, drop_first=True)
Теперь посмотрим:
dt.head()
Хорошо, теперь перейдем к модели.
Следующая часть подгоняет модель случайного леса к данным:
X_train, X_test, y_train, y_test = train_test_split(dt.drop('target', 1), dt['target'], test_size = .2, random_state=10) #split the data
model = RandomForestClassifier(max_depth=5)
model.fit(X_train, y_train)
Мы можем построить последовательное дерево решений, чтобы увидеть, что оно делает:
estimator = model.estimators_[1]
feature_names = [i for i in X_train.columns]
y_train_str = y_train.astype('str')
y_train_str[y_train_str == '0'] = 'no disease'
y_train_str[y_train_str == '1'] = 'disease'
y_train_str = y_train_str.values
#код из https://towardsdatascience.com/how-to-visualize-a-decision-tree-from-a-random-forest-in-python-using-scikit-learn-38ad2d75f21c
export_graphviz(estimator, out_file='tree.dot',
feature_names = feature_names,
class_names = y_train_str,
rounded = True, proportion = True,
label='root',
precision = 2, filled = True)
from subprocess import call
call(['dot', '-Tpng', 'tree.dot', '-o', 'tree.png', '-Gdpi=600'])
from IPython.display import Image
Image(filename = 'tree.png')
Это дает нам инструмент для объяснения. Однако сейчас не возможно получить представление о наиболее важных функциях. Поэтому мы вернемся к ним позже. Давайте оценим модель:
y_predict = model.predict(X_test)
y_pred_quant = model.predict_proba(X_test)[:, 1]
y_pred_bin = model.predict(X_test)
Оценим соответствие с помощью матрицы ошибок:
confusion_matrix = confusion_matrix(y_test, y_pred_bin)
confusion_matrix
Диагностические тесты часто используются с sensitivity(чувствительностью) и specificity(специфичностью) в качестве основных показателей. Чувствительность и специфичность определяются как:
Давайте посмотрим, что дает эта модель:
total=sum(sum(confusion_matrix))
sensitivity = confusion_matrix[0,0]/(confusion_matrix[0,0]+confusion_matrix[1,0])
print('Sensitivity : ', sensitivity )
specificity = confusion_matrix[1,1]/(confusion_matrix[1,1]+confusion_matrix[0,1])
print('Specificity : ', specificity)
Выглядит разумным. Давайте также проверим с помощью Receiver Operator Curve (ROC),
fpr, tpr, thresholds = roc_curve(y_test, y_pred_quant)
fig, ax = plt.subplots()
ax.plot(fpr, tpr)
ax.plot([0, 1], [0, 1], transform=ax.transAxes, ls="--", c=".3")
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])
plt.rcParams['font.size'] = 12
plt.title('ROC curve for diabetes classifier')
plt.xlabel('False Positive Rate (1 - Specificity)')
plt.ylabel('True Positive Rate (Sensitivity)')
plt.grid(True)
Другой распространенный показатель - Площадь под кривой или AUC. Это удобный способ зафиксировать производительность модели одним числом, хотя и не без определенных проблем. Как показывает практика, AUC можно классифицировать следующим образом:
Посмотрим, что нам дает вышеупомянутая ROC:
auc(fpr, tpr)
Хорошо, значит, все работает хорошо.
Теперь давайте посмотрим, что дает нам модель с помощью инструментов объяснимости машинного обучения.
Permutation importance(Важность перестановки) это первый инструмент для понимания модели машинного обучения, который включает перетасовку отдельных переменных в данных проверки (после того, как модель была подогнана), и наблюдение за влиянием на точность. Узнайте больше here.
Давайте взглянем:
perm = PermutationImportance(model, random_state=1).fit(X_test, y_test)
eli5.show_weights(perm, feature_names = X_test.columns.tolist())
Итак, похоже, что наиболее важным фактором с точки зрения перестановки является талесемия в результате «обратимого дефекта». Высокая важность типа «достигнутая максимальная частота сердечных сокращений» имеет смысл, поскольку это непосредственное субъективное состояние пациента во время обследования (в отличие, скажем, от возраста, который является гораздо более общим фактором).
Давайте внимательнее посмотрим на количество основных сосудов, используя Partial Dependence Plot(График частичной зависимости) (here). Эти графики изменяют одну переменную в одной строке в диапазоне значений и смотрят, как она влияет на результат. Он делает это для нескольких строк и отображает средний эффект. Давайте посмотрим на переменную 'num_major_vessels', которая была в верхней части списка важности перестановок:
base_features = dt.columns.values.tolist()
base_features.remove('target')
feat_name = 'num_major_vessels'
pdp_dist = pdp.pdp_isolate(model=model, dataset=X_test, model_features=base_features, feature=feat_name)
pdp.pdp_plot(pdp_dist, feat_name)
plt.show()
Итак, мы можем видеть, что по мере увеличения числа крупных кровеносных сосудов вероятность сердечных заболеваний уменьшается. Это имеет смысл, поскольку означает, что больше крови может попасть в сердце.
А что насчет возраста?
feat_name = 'age'
pdp_dist = pdp.pdp_isolate(model=model, dataset=X_test, model_features=base_features, feature=feat_name)
pdp.pdp_plot(pdp_dist, feat_name)
plt.show()
Это немного странно. Чем старше возраст, тем меньше вероятность сердечных заболеваний? Хотя синие доверительные области показывают, что это может быть неправдой (красная базовая линия находится в синей зоне).
А насчет 'st_depression':
feat_name = 'st_depression'
pdp_dist = pdp.pdp_isolate(model=model, dataset=X_test, model_features=base_features, feature=feat_name)
pdp.pdp_plot(pdp_dist, feat_name)
plt.show()
Интересно, что эта переменная также показывает уменьшение вероятности, чем выше она идет. Что именно это? Поиск привел меня к следующему описанию Энтони Л. Комарова, доктора медицины, специалиста по внутренним болезням 5 .... "Электрокардиограмма (ЭКГ) измеряет электрическую активность сердца. Волны, которые появляются на ней, обозначаются буквами P, QRS и T. Каждая из них соответствует разным частям сердцебиения. ST сегмент представляет электрическую активность сердца сразу после сокращения правого и левого желудочков, перекачивая кровь к легким и остальным частям тела. После этого большого усилия мышечные клетки желудочков расслабляются и готовятся к следующему сокращению. В течение этого периода, электричество проходит мало или отсутствует, поэтому сегмент ST находится на одном уровне с базовой линией или иногда немного выше нее. Чем быстрее сердце бьется во время ЭКГ, тем короче становятся все волны. Форма и направление ST сегмент гораздо важнее его длины. Вверх или вниз смещение палаты может означать снижение притока крови к сердцу по разным причинам, включая сердечный приступ, спазмы в одной или нескольких коронарных артериях (стенокардия Принцметала), инфекция внутренней оболочки сердца (перикардит) или самой сердечной мышцы (миокардит), избыток калия в кровотоке, нарушение сердечного ритма или сгусток крови в легких (тромбоэмболия легочной артерии)».
Таким образом, эта переменная, которая описывается как «депрессия ST, вызванная упражнениями по сравнению с отдыхом», кажется, предполагает, что чем выше значение, тем выше вероятность сердечного заболевания, но график выше показывает обратное. Возможно, важна не только величина депрессии, но и взаимодействие с типом наклона? Давайте проверим с помощью 2D PDP:
inter1 = pdp.pdp_interact(model=model, dataset=X_test, model_features=base_features, features=['st_slope_upsloping', 'st_depression'])
pdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=['st_slope_upsloping', 'st_depression'], plot_type='contour')
plt.show()
inter1 = pdp.pdp_interact(model=model, dataset=X_test, model_features=base_features, features=['st_slope_flat', 'st_depression'])
pdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=['st_slope_flat', 'st_depression'], plot_type='contour')
plt.show()
Похоже, что в обоих случаях низкая депрессия - это плохо.
Посмотрим, что нам говорят значения SHAP. Они работают, показывая влияние значений каждой переменной в одной строке по сравнению с их базовыми значениями. (here).
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)
shap.summary_plot(shap_values[1], X_test, plot_type="bar")
Количество крупных сосудов свверху. Воспользуемся сводным графиком значений SHAP:
shap.summary_plot(shap_values[1], X_test)
Количество основных судов довольно ясно, и это говорит о том, что низкие значения - это плохо (синий справа). Разделение талассемии на «обратимый дефект» очень четкое (да = красный = хорошо, нет = синий = плохо).
Вы можете увидеть четкое разделение многих других переменных. Стенокардия, вызванная физической нагрузкой, имеет четкое разделение, хотя и не так, как ожидалось, поскольку «нет» (синий) увеличивает вероятность. Другой очевидный - st_slope. Похоже, когда он плоский, это плохой знак (красный справа).
Также странно то, что у мужчин (красный) в этой модели снижен шанс сердечных заболеваний. Почему это? Знание предметной области говорит нам, что у мужчин больше шансов.
Затем давайте выберем отдельных пациентов и посмотрим, как различные переменные влияют на их результаты.
def heart_disease_risk_factors(model, patient):
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(patient)
shap.initjs()
return shap.force_plot(explainer.expected_value[1], shap_values[1], patient)
data_for_prediction = X_test.iloc[1,:].astype(float)
heart_disease_risk_factors(model, data_for_prediction)
Для этого человека их прогноз составляет 36% (по сравнению с исходным уровнем 58,4%). Многие вещи работают в их пользу, в том числе наличие крупного сосуда, обратимого дефекта талассемии и отсутствия плоского наклона.
Давай проверим еще:
data_for_prediction = X_test.iloc[3,:].astype(float)
heart_disease_risk_factors(model, data_for_prediction)
Для этого человека их прогноз составляет 70% (по сравнению с исходным уровнем 58,4%). Не работают в их пользу такие вещи, как отсутствие крупных сосудов, плоский наклон и не обратимый дефект талассемии.
Мы также можем построить так называемые «графики вклада зависимостей SHAP» (here), которые довольно очевидны в контекст значений SHAP:
ax2 = fig.add_subplot(224)
shap.dependence_plot('num_major_vessels', shap_values[1], X_test, interaction_index="st_depression")
Вы можете увидеть резкий эффект на количество крупных сосудов, но, похоже, цвет (st_depression) не дает многого.
Заключительный график - один из самых эффективных. Он показывает прогнозы и факторы влияния для многих (в данном случае 50) пациентов, вместе взятых. Наведите указатель мыши, чтобы увидеть, почему каждый человек оказался красным (предсказание болезни) или синим (предсказание отсутствия болезни):
shap_values = explainer.shap_values(X_train.iloc[:50])
shap.force_plot(explainer.expected_value[1], shap_values[1], X_test.iloc[:50])
Этот набор данных старый и маленький по сегодняшним меркам. Однако это позволило создать простую модель, а затем использовать различные инструменты и методы объяснения машинного обучения, чтобы заглянуть внутрь. Вначале было построено предпложение, используя общие знания в области, что такие факторы, как холестерин и возраст, будут основными факторами в модели. Этот же набор данных не показал этого. Вместо этого преобладали основные факторы и аспекты результатов ЭКГ.
Я подозреваю, что такой подход будет приобретать все большее значение, поскольку с каждый годом все больше растут вычислительные мощности, а как следсвие и машинное обучение играет все более важную роль в нашей жизни, в том числе и в здравоохранении.